function [B,Elist,Eulist] = demons_optimizer(M,F,AFT,T,sizes,iscomp,options)

%% Check/set input options
defaultoptions=struct('sigma_fluid',8, 'sigma_diff',0,'alpha',1.5,'interpolation','linear','maxiter',500);
if(~exist('options','var')),
    options=defaultoptions;
else
    tags = fieldnames(defaultoptions);
    for i=1:length(tags)
        if(~isfield(options,tags{i})),  options.(tags{i})=defaultoptions.(tags{i}); end
    end
    if(length(tags)~=length(fieldnames(options))),
        warning('BsplineRegistrationGradient:unknownoption','unknown options found');
    end
end

% Set all parameters
alpha=options.alpha;
maxiter = options.maxiter;
sigma_diff=options.sigma_diff;
sigma_fluid=options.sigma_fluid;
if(strcmpi(options.interpolation(1),'l')),
    interpolation_mode=0;
else
    interpolation_mode=2;
end

%%

S = reshape(T, sizes);
sx = S(:,:,1); sy = S(:,:,2);
sxpre = sx; sypre = sy;

ux = zeros(size(sx,1),size(sx,2));
uy = zeros(size(sy,1),size(sy,2));

% Calculate current moving image
[X,Y] = ndgrid(0:(sizes(1)-1),0:(sizes(2)-1));
Xd=X-(sizes(1)/2); Yd=Y-(sizes(2)/2);

AFx = ((sizes(1)/2) + AFT(1,1) * Xd + AFT(1,2) *Yd + AFT(1,3) * 1)-X;
AFy = ((sizes(2)/2) + AFT(2,1) * Xd + AFT(2,2) *Yd + AFT(2,3) * 1)-Y;
[TBx, TBy] = VFcomposition(AFx, AFy, sx+ux, sy+uy, size(sx,2), size(sx,1));

MAF = movepixels(M, AFx, AFy, [], interpolation_mode);

%%
iscomp = 1;
if iscomp == 1
    Ms = movepixels(M,TBx,TBy,[],interpolation_mode);
else
    Ms = movepixels(MAF,sx+ux,sy+uy,[],interpolation_mode);
end
%%
% figure, imshow(Ms);
% Calculate difference image
Idiff=Ms-F;

% Calculate error after current transformation
Em =Idiff.^2+alpha^2*Idiff.^2.*(ux.^2+uy.^2);
Emu = alpha^2*Idiff.^2.*(ux.^2+uy.^2);
E = sum(Em(:));
Eu = sum(Emu(:));

Eold = E;
Ediff = -100;

disp(['original energy = ', num2str(E)]);
Elist(1) = E;
Emlist(1) = Eu;

iteration = 0;
%%
while 1
    
    %     Ms = movepixels(MAF,sx+ux,sy+uy,[],interpolation_mode);
    %     % Calculate difference image
    %     Idiff=Ms-F;
    %
    %     % Calculate error after current transformation
    %     Em =Idiff.^2+alpha^2*Idiff.^2.*(ux.^2+uy.^2);
    
    %     E = sum(Em(:));
    
    % Calculate gradients in both images
    [My,Mx] = gradient(Ms);
    [Fy,Fx] = gradient(F);
    
    % Calculate gradient (update field) in x and y direction
    Ma=1./(eps+(Mx.^2+My.^2)+alpha^2*Idiff.^2);
    Fa=1./(eps+(Fx.^2+Fy.^2)+alpha^2*Idiff.^2);
    ux=-(Idiff.*(Fx.*Fa)+Idiff.*(Mx.*Ma));
    uy=-(Idiff.*(Fy.*Fa)+Idiff.*(My.*Ma));
    
    % Smooth gradient incase of fluid-like regularization
    if(sigma_fluid>0)
        ux=imgaussian(ux,sigma_fluid);
        uy=imgaussian(uy,sigma_fluid);
    end
    
    cx = sx + ux;
    cy = sy + uy;
    
    % Smooth gradient incase of diffusion-like regularization
    if(sigma_diff>0)
        % Smooth transformation field
        sx=imgaussian(cx,sigma_diff);
        sy=imgaussian(cy,sigma_diff);
    else
        sx=cx; sy=cy;
    end
    [TBx, TBy] = VFcomposition(AFx, AFy, sx+ux, sy+uy, size(sx,2), size(sx,1));
    
    if iscomp == 1
        Ms = movepixels(M,TBx,TBy,[],interpolation_mode);
    else
        Ms = movepixels(MAF,sx+ux,sy+uy,[],interpolation_mode);
    end
    
    %     figure, imshow(Ms);
    % Calculate difference image
    Idiff=Ms-F;
    % Calculate error after current transformation
    Em =Idiff.^2+alpha^2*Idiff.^2.*(ux.^2+uy.^2);
    Emu = alpha^2*Idiff.^2.*(ux.^2+uy.^2);
    E = sum(Em(:));
    Eu = sum(Emu(:));
    Ediff = E - Eold;
    Eold = E;
    %     if mod(iteration,25)==0
    %         close all
    %     end
    iteration = iteration+1;
    Elist(iteration+1) = E;
    Eulist(iteration+1) = Eu;
    
    if iteration>=maxiter
        break;
    else if Ediff>=0
            sx = sxpre; sy = sypre;
            Elist(iteration+1) = [];
            Eulist(iteration+1) = [];
            break;
        end
    end
    
    disp(['iteration = ', num2str(iteration), ', energy = ', ...
        num2str(E), ', energy diff = ', num2str(Ediff)] );
    sxpre = sx; sypre = sy;
end
% Both x and y gradient to one output variable

% Egrad=zeros(size(T));
% Egrad(:,:,1)=ux;
% Egrad(:,:,2)=uy;
% Egrad=Egrad(:);
% X = 1:size(Elist,2);
% plot(X, Elist,'-r*' );
Bm(:,:,1) = sx; Bm(:,:,2) = sy;
B = Bm(:);
end
